﻿using System;
using System.Text;
using System.Drawing;
using System.Buffers;
using System.Collections;
using System.Collections.Generic;
using System.Runtime.InteropServices;

public static partial class NativeAOT
{
    [UnmanagedCallersOnly(EntryPoint = "lluu")]
    public static unsafe int lluu(IntPtr a_ptr, int n, IntPtr l_ptr, IntPtr u_ptr)
    {
        double* a = (double*)a_ptr.ToPointer();
        double* l = (double*)l_ptr.ToPointer();
        double* u = (double*)u_ptr.ToPointer();

        return lluu(a, n, l, u);
    }

    /// <summary>
    /// 矩阵的三角分解
    /// </summary>
    /// <param name="a">a[n][n]存放n阶矩阵A。返回时存放Q矩阵。</param>
    /// <param name="n"></param>
    /// <param name="l">l[n][n]返回时存放下三角矩阵L。</param>
    /// <param name="u">u[n][n]返回时存放上三角矩阵U。</param>
    /// <returns>若为0，则表示失败；若不为0，则表示正常。</returns>
    public static unsafe int lluu(double* a, int n, double* l, double* u)
    {
        int i, j, k, w, v, ll;

        for (k = 0; k <= n - 2; k++)
        {
            ll = k * n + k;
            //右下角子阵中a[k][k]=0
            //if (Math.Abs(a[ll]) <= float.Epsilon)
            if (EQ0(a[ll]))
            {
                return 0;
            }
            for (i = k + 1; i <= n - 1; i++)
            {
                w = i * n + k; a[w] = a[w] / a[ll];
            }
            for (i = k + 1; i <= n - 1; i++)
            {
                w = i * n + k;
                for (j = k + 1; j <= n - 1; j++)
                {
                    v = i * n + j;
                    a[v] = a[v] - a[w] * a[k * n + j];
                }
            }
        }
        for (i = 0; i <= n - 1; i++)
        {
            // L-I
            for (j = 0; j < i; j++)
            {
                w = i * n + j;
                l[w] = a[w];
                u[w] = 0.0;
            }
            w = i * n + i;
            l[w] = 1.0;
            u[w] = a[w];
            // U
            for (j = i + 1; j <= n - 1; j++)
            {
                w = i * n + j;
                l[w] = 0.0;
                u[w] = a[w];
            }
        }
        return 1;
    }

#if __DRIVE_CODE__
      int main()
      { 
          int i,j;
          double l[4][4],u[4][4];
          double q[4][4],a[4][4]={ {2.0,4.0,4.0,2.0},
                                {3.0,3.0,12.0,6.0},
                                {2.0,4.0,-1.0,2.0},
                                {4.0,2.0,1.0,1.0}};
          for (i=0; i<=3; i++)
              for (j=0; j<=3; j++)  q[i][j]=a[i][j];
          i=lluu(&q[0][0],4,&l[0][0],&u[0][0]);
          if (i==0) return 0;
          cout <<"MAT Q=L+U-I:\n";
          for (i=0; i<=3; i++)
          {
              for (j=0; j<=3; j++)  cout <<q[i][j] <<"    ";
              cout <<endl;
          }
          cout <<"MAT L:\n";
          for (i=0; i<=3; i++)
          {
              for (j=0; j<=3; j++)  cout <<l[i][j] <<"    ";
              cout <<endl;
          }
          cout <<"MAT U:\n";
          for (i=0; i<=3; i++)
          { 
              for (j=0; j<=3; j++)  cout <<u[i][j] <<"    ";
              cout <<endl;
          }
          return 0;
      }
#endif
}